Validation of SNP markers for thermotolerance adaptation in Ovis aries adapted to different climatic regions using KASP-PCR technique

A study on 51 SNPs belonging to 29 genes related to heat stress was carried out in 720 sheep from 17 different breeds adapted to different climates from Hungary, Bosnia and Herzegovina, Morocco and Romania, using Kompetitive Allele-Specific Polymerase Chain Reaction. Genotype frequency and the Hardy–Weinberg equilibrium were calculated, followed by a clustering using the Principal Component Analysis. We analyzed the polymorphisms in the following genes analyzed: HSPA12A, HSP90AA1, IL33, DIO2, BTNL2, CSN2, ABCG1, CSN1S1, GHR, HSPA8, STAT3, and HCRT. We emphasized on HSPA12A and HSPA8 genes as they were successfully genotyped in all studied flocks in which genotype frequency patterns were identified. Contrary to previous findings, the A allele for HSPA8 SNP was not observed in the heat tolerant breeds, being found exclusively in cold-tolerant breeds. The principal component analysis could not clearly differentiate the breeds, while plot concentration was slightly varied among the three groups, with HSP90AA1 and IL33 SNPs’ loading values significantly contributing to PC1 and PC2. We confirmed previous works that the HSPA12A, HSPA8, HSP90AA1 and IL33 SNPs are potential candidate markers for thermotolerance adaptation in sheep. This research contributes to the genetic variability of SNPs for thermotolerance adaptability in sheep.

The heterozygote TC for rs161504783-HSPA12A was dominant, except for Hungarian Racka, Transylvanian Merino, Hungarian Merino, Botosani Karakul and Sardi. The T and C allele frequencies were almost equally frequent in most breeds, except in Sardi (T allele = 0.107 and C alleles = 0.828) and Botosani Karakul (T allele = 0.735 and C alleles = 0.265).
The heterozygote GA for SNP rs588145625-HSPA8 was absent in the heat tolerance breeds, except the Transylvanian Merino, and in breeds with high cold tolerance (Hungarian Racka, Babolna Tetra, Hungarian Tsigai, Romanian Racka, Pramenka and Turcana). The homozygote GG was present in all heat-tolerant breeds and some cold-tolerant breeds (Suffolk, Ile de France and Hungarian Merino). The G allele was dominant in all breeds, ranging from 0.760 to 1, and the A allele has only appeared in cold-tolerant breeds and Transylvanian Merino with frequency ranging from 0 to 0.308.
No patterns have been found for allelic and genotypic frequencies for rs588498137-STAT3 and rs602521720-HCRT SNPs.
The GG genotype and G allele were dominant in all breeds for SNP rs588498137-STAT3, with G allelic frequency varying from 0.750 to 1. Similarly, the CC genotype and C allele for SNP rs602521720-HCRT was dominant in all breeds, with C allelic frequency varying from 0.546 to 1.
The Hardy-Weinberg equilibrium (HWE) was investigated with Chi-square (x 2 ) test using allelic frequencies, observed and expected genotypes, and P value of the polymorphic genes are summarized in Table S2. Botosani Karakul breed from Romania was found to be the breed with most SNPs deviated from HWE (P ≤ 0.05); rs416259751-IL33, rs411181557-DIO2, rs414917134-BTNL2 and rs420611298-ABCG1.
Genetic diversity and interrelationship between SNPs. SNPs data were used to perform Principal Component Analysis (PCA) to highlight differences between either the breeds (Fig. 1a) or the climatic characteristics (Fig. 1b).

Discussion
The genetic background of thermotolerance adaptation in sheep represents a great interest for current research, especially when faced with the imminent increase of the earth's temperature because of climate change. Different breeds from different regions were used in this study, with the aim to see the different genetic backgrounds of thermotolerance adaptation according to geographical origins. The application of KASP-PCR assay in this study did not show a high assay success rate as in previous studies carried out on goats 28 and dairy cattle 29 . As much as 62.74% of SNPs were successfully genotyped in this study. These relatively high (37.26%) unsuccessful genotypes could be the consequence of sample damage during transportation from the laboratory to the subcontracting laboratory. Furthermore, quality control results in lower data for the analysis.
Only four SNPs were successfully genotyped in all investigated populations, and the focus was given to SNPs HSPA12A and HSPA8 as they are the only ones that showed allelic and genotypic frequency patterns with the climatic characteristics. HSPs are a wide family of chaperone proteins categorized according to their molecular size and amino acid sequence similarity. The HSPA12A and HSPA8 are members of the HSP70 family with a molecular weight of 70 kDa 30 . It is the biggest, most numerous and most conserved protein family throughout evolution, as well as the most extensively studied protein family across a wide range of species [31][32][33] . It is found on chromosome 15 and is made up of nine exons in sheep. This gene is widely studied for thermal adaptability. It controls cellular survival to heat stress with a highly dynamic nature, and these proteins are responsible for maintaining the organism's equilibrium and acclimating to heat stress 34,35 .
HSPA12A affects aspects like cellular senescence and how cells respond to heat stress 36 . It was found to be more active in ruminants during the summer, which helps them adapt to a hostile environment. This gene was found to be highly associated with heat tolerance in indicine cattle (Bos taurus indicus) and also has an important role in the water holding capacity in beef breeds, which might also act as a mechanism in surviving heat stress through the increased adrenaline and pH change in muscles 37,38 . According to the genotype and allele frequency results, the heterozygote TC for rs161504783-HSPA12A was dominant in most of the breeds, except Hungarian Racka, Hungarian Merino, Transylvanian Merino, Botosani Karakul and Sardi. Béni Guil, D'Man and Timahdite, breeds from a hot region, were not in HWE for HSPA12A.
In HSP70, exon 1 is non-coding, while the following eight exons combine to form the HSPA8 protein, with 650 amino acids of 71 kDa. HSPA8 aids in the day-to-day cell functions of protein folding and unfolding, polypeptide aggregation prevention, disassembly of large protein complexes and protein translocation across cellular compartments 35,39 . Given that increasing HSPA8 levels have been discovered to be positively correlated with heat tolerance, this gene has been employed as a candidate for heat resistance in many livestock species [40][41][42] .
In this study, the homozygote GG was carried by all heat-tolerant breeds and some cold-tolerant breeds (Suffolk, Ile de France and Romanian Tsigai), while the heterozygote GA for SNP rs588145625-HSPA8 was only found in cold-tolerant breeds (Hungarian Racka, Babolna Tetra, Hungarian Tsigai, Romanian Racka, Transylvanian Merino, Pramenka and Turcana). The A allele was absent in heat tolerant breeds, which is opposed to the observation made by 25 in a study of gene expression in Indian sheep, which discovered that animals with A allele have better hot climate adaptability, compared to animals with the G allele. The AA genotype is more adaptable to a hot environment and has lower HSPA8 gene expression, than animals with the AG genotype. Similarly, a study by 41 showed that the GG genotype has the least ability to survive heat stress in Awassi and Arabi sheep.
In this study, Botosani Karakul has shown to be the breed with the largest number of SNPs deviating from HWE. Lower heterozygosity was observed in 10 SNPs, and 4 of them deviated from HWE. One possible reason for this deviation of HWE was because Botosani Karakul is one of the Romanian crossbreeds that has undergone extensive genetic mixing from its German and Austria lines, and other Romanian breeds since its introduction from Russia at the beginning of the nineteenth century 43,44 .
PCA was unable to differentiate each breed in a clear and distinct manner; however, a relative clustering was observed based on the three different categories of climatic regions (Fig. 1b). This could be attributed to the fact that heat-tolerant breeds that are kept in the EU have become acclimated to the subtropical environment, in addition to the fact that genetic ad-mixture has occurred as a result of the widespread use of reproductive technologies, which has led to a less distinct genetic divergence between breeds, with previous studies reporting high levels of ad-mixture between different sheep breeds reared in Hungary and Romania, with the aim to increase the productivity [45][46][47][48][49] . We acknowledge that the cold breeds used in this study were not from year-around cold regions (e.g. Iceland, Finland, Norway), which was one of our limitations due to our inability to obtain samples from these regions. However, we strongly believe that the cold tolerance breeds from temperate climates utilized in this study are sufficiently contrasting with the Moroccan heat tolerant breeds.
From the loading biplot, HSP90AA1 and IL33 SNPs significantly contributed to PC1 and PC2. HSP90AA1 has been confirmed to be associated with thermal stress susceptibility of sheep in previous studies 50,51 , while IL33 was found to affect the sheep immunity and resistance to gastrointestinal intestinal nematode infection 52,53 , which is also associated because heat stress promotes to immune suppression and increases animal vulnerability to illnesses 54 . Our PCA findings validate these two SNPs as potential candidates for heat adaptability across the sheep breeds investigated in the current study, although more research is needed in order to clarify the relationship between the investigated SNPs and heat resistance.

Conclusion
Based on our 17 SNP polymorphism analyses performed on 601 animals, we found that the KASP-PCR method represents a feasible method for investigating polymorphisms in different sheep breeds. Furthermore, based on allele and genotype frequency, we validated that HSPA12A and HSPA8 SNPs are potential candidate markers for thermotolerance adaptation in sheep, whereas principal component analysis confirmed that HSP90AA1 and www.nature.com/scientificreports/ IL33 SNPs were the primary potential candidates. The results contribute to an increase in knowledge regarding the genetic variability of SNPs for thermotolerance adaptation in sheep. However, more studies are needed in order to clarify the relationship between the studied SNPs and heat resistance in sheep.

Method
Genomic DNA extraction. Samples were collected from 720 sheep belonging to17 breeds adapted to different climatic conditions, originating from 4 countries ( Table 1). The majority of the breeds studied are indigenous to the country of origin, with some exotic breeds being included to determine if their acclimatization has contributed to adaptation. The breed characteristic (hot or cold tolerant) was determined based on the origin and development history of each breed in that particular country, as well as the environmental conditions under which the sample was collected (  Genotyping and quality control. The bi-allelic discrimination of the selected 51 SNPs was performed using Kompetitive Allele Specific PCR (KASPTM, LGC Genomics, Teddington, Middlesex, UK). SNP Viewer software version 1.99 (Hoddesdon, UK) was used to visualize the results. All genotype data were exported for statistical analysis. Only SNPs that appeared in at least 50% of the breeds were considered. Data quality control of genotyped data included discarding animals with a call rate of less than 50% and the SNPs with call rates < 50%. This led to discrepancies in either the number of animals per breed or the number of SNPs per animal.
Data analysis. The raw allele calls obtained from LGC Genomics were analyzed using LGC Genomics' KlusterCaller program. Gene diversity, allele and genotype frequencies, and their accordance with or deviation from the Hardy-Weinberg equilibrium were determined by POPGENE software version 1.31 60 The Principal Component Analysis was done using FactoMineR 61 and ggplot2 62 packages from the R Program 63 to visualize the genetic divergences between sheep breeds that were divided according to their climatic characteristics; cold-tolerant breeds (Babolna Tetra, Hungarian Merino, Hungarian Racka, Hungarian Tsigai, Ile de France, Pramenka, Romanian Racka, Suffolk and Turcana), heat tolerance breeds originated from Morocco (Béni Guil, D'Man, Timahdite and Sardi), and heat tolerant breeds reared in Europe (Hungarian Awassi, Botosani Karakul, Transylvanian Merino, and Romanian Tsigai).
Ethical approval. The authors confirm that the experiment complied with the European Union's Directive on Animal Experimentation (Directive 2010/63/EU) and ARRIVE guidelines. All animals in the experiment underwent standard procedures without experiencing any harm or discomfort, and all procedures were carried out in compliance with applicable guidelines and regulations. The study was approved by the Scientific and Ethics Committee of Centre for Agricultural Genomics and Biotechnology, University of Debrecen (Ethics statement No. 07).

Data availability
The 10 SNPs genotype datasets generated and analyzed in this study are only partially available at the EVA under accession number ERZ6760182. The remaining 7 SNPs genotype data is available in the supplementary material file (Table S6)